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^ ■ Abstract 

We propose a new model of turbulence for use in large-eddy simulations (LES). The turbulent force, 

represented here by the turbulent Lamb vector, is divided in two contributions. The contribution including 

only subfilter fields is deterministically modeled through a classical eddy- viscosity. The other contribution 

including both filtered and subfilter scales is dynamically computed as solution of a generalized (stochastic) 

Langevin equation. This equation is derived using Rapid Distortion Theory (RDT) applied to the subfilter 

scales. The general friction operator therefore includes both advection and stretching by the resolved scale. 

The stochastic noise is derived as the sum of a contribution from the energy cascade and a contribution 
CZ2 ■ . ... 

from the pressure. The LES model is thus made of an equation for the resolved scale, including the 

turbulent force, and a generalized Langevin equation integrated on a twice-finer grid. We compare the 

full model with several approximations. In the first one, the friction operator of the Langevin equation 

• is simply replaced by an empirical constant, of the order of the resolved scale correlation time. In the 

second approximation, the integration is replaced by a condition of instantaneous adjustment to the 

stochastic force. In this approximation, our model becomes equivalent to the velocity- estimation model 

y— I , of Domaradzki et al. [ZHEHHj- In the isotropic, homogeneous situations we study, both approximations 

provide satisfactory results, at a reduced computational cost. The model is finally validated by comparison 

to DNS and is tested against classical LES models for isotropic homogeneous turbulence, based on eddy 

viscosity. We show that even in this situation, where no walls are present, our inclusion of backscatter 

, through the Langevin equation results in a better description of the flow. 

o ■ 
\o ■ 

O '. 1 Introduction 

The rationale for Large Eddy Simulation is often rooted in our inability to handle all the degrees of 
freedom of a large Reynolds number turbulent flow. Given that the smaller scales monopolize most of 
the computing resources, it is tempting to cut through the number of degrees of freedom via an ad hoc 
. small-scale decimation. The price to pay is of course a need for parameterization (the so-called SubGrid 

Scale models, or SGS), to make up for the energy transfer of the ghost scales. We refer the reader to 
01321 HH EH for recent reviews about modern parameterization strategies in LES. We may loosely divide 
the SGS models in two classes, according to their philosophy: the "deterministic" models based on eddy- 
viscosities, and the "stochastic" models, based on synthetic fields. In eddy-viscosity methods, the action 
of small scales is parameterized via a few deterministic numbers, linked with the various components of 
the subgrid-scale stress tensor. These models seek to reproduce the intensification of energy transport 
due to the action of scales widely separated from the considered one. However, they fail to reproduce 
backward energy transfer (backscatter) from small to large scale, created by elongated triads in the 
spectral spaceEO]. This effect has been shown to induce a stochastic character in the LES |2Jj]. In ideal 
situations, where the turbulence is isotropic, homogeneous and far from wall, this backscatter is usually 
viewed as secondary, and eddy-viscosity based models are generally satisfactory. However, in more realistic 
situations, including turbulence near walls or in boundary layers, the energy backscatter has proven to 
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reproduce a realistic backscatter in some situations ( [151 I23[ ). The need for backscatter modeling also 
leads to the development of "stochastic" strategies, where the discarded small-scale motions are replaced 
by a set of random numbers, mimicking cither a random force or synthetic velocity fields [71 1331 Il7| . In 
many ways, these strategies resemble the strategies used to describe the dynamics of a heavy particle 
coupled to a thermal bath involving many degrees of freedom (the so-called Brownian motion). The 
decimation is here performed by substituting in place of the bath a deterministic friction and a stochastic 
force, the two terms being linked through the dissipation theorem. The initial problem is then completely 
described through the so-called Langevin equation. 

Turbulence is typically an out-of-equilibrium system, and there is probably no hope that such a simple 
description will ever be possible (would it be only because no fluctuation-dissipation theorem prevails in 
turbulence!). However, we would like to use this analogy to motivate a new strategy for LES modeling: 
replace the actual dynamics of the decimated degrees of freedom by a suitable noise, via a Langevin 
equation. Although this strategy may seem close to recent models based on synthetic fields, we would like 
to point out an important philosophical difference: rather than trying to estimate the actual small-scale 
dynamics, we aim at trying to estimate a plausible small-scale dynamics. We believe there is no unique 
solution for this last option. In the sequel, we present one solution based upon Rapid Distortion Theory 
[17] . There may exist actually more efficient models, based e.g. upon information theory |lfi[ . 

To be more specific, consider a turbulent flow, with velocity field Ui(x,t) and introduce a filtering 
procedure so as to separate it into a resolved field Ui = u[ and a subfilter field u\ = Ui — J7,. The 
resolved field obeys a dynamical equation obtained by filtering of the Navier-Stokes equation, which may 
conveniently be written as |39j : 



Here, P is the resolved pressure and v is the viscosity and 1 is a turbulent force . In idealistic situations 
including a spectral gap between resolved and subfilter scale, the vector 1 is zero, and one can rigorously 
show that the contribution of the remaining term is of " diffusive" type (providing certain symmetries which 
exclude first order behavior such as the AKA effect IT ). Even as one departs from this idealistic situation, 
experimental [27J and numerical study ^Hj show that this term correlates strongly to the resolved velocity 
gradient, thereby allowing a deterministic treatment through an eddy- viscosity of appropriate shape. In 
the same time, the force vector 1 becomes increasingly non-negligible (it can even become dominant in 2D 
situations see ^JJ). It is responsible for backscatter type of behaviour and needs to be modeled through 
novel "non-diffusive" and "non-deterministic" strategies. In the sequel, we will therefore focus onto the 
modeling of the 1 term, via a generalized Langevin equation <9 t l = ^41 + £ where A is a generalized evolution 
operator, and £ is a noise. 

However, a clear difficulty associated with this strategy is the lack of theoretical guide (equivalent 
of the statistical mechanics in Brownian motion) to help us devising the "best" friction, and the "best" 
stochastic force. For the time being, we then choose to pin down our model as best as possible to the real 
dynamics of Navier-Stokes by trying to derive it from the original dynamical equations, rather than from 
pure empirical or dimensional considerations. For this, we reformulate the RDT-Langevin model of Laval 
et al. 17 in a way suitable for LES. There are of course limitations to this approach, pertaining the need 
for both a simple enough model, and for tractable analytical computations. We try to formulate them 
as honestly as possible by pointing out the approximation we make at the various stage of the derivation 
of the model. This is done in Section The LES-Langevin model is then implemented in the case 
of periodic three-dimensional turbulence and comparison with other existing LES models is provided in 
Section Our conclusions follow in Section ^] 

2 The LES-Langevin model of turbulence 

2.1 Derivation of the Langevin equation 

Our derivation is based on the stochastic RDT model developed by Laval et al. ^UEIE]- This model is 
based on the observation that subfilter-scales are mostly slaved to resolved scales via linear processes akin 
to rapid distortion. This property is substantiated by various numerical simulations, and is linked with 
the prominence of non-local interactions at subfilter-scale 117} . Using incompressibility, the small-scale 
dynamics in this model can be written as: 



dtU + (U • V)U + 1 + (u' • V)u' 
1 = (U • V)u' + (u' • V)U. 



VP + j/AU 



(1) 



d t u' 



¥ - Vp' + V(v + W()Vu' - f . 



(2) 
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subfilter scales and f is a forcing stemming from the energy cascade. The latter can be shown to be 
dominated by resolved scales non-linearities through /, = dj \U%Uj — UiUj). Finally, we may use the 
observation that subfilter scales vary over fast time scale with respect to resolved scales to write: 

d t l m (U • V)<9tu' + (d t u' ■ V)U, 

« -{(U-V)(l' + f)) + [(I' + f)-V]U} 
- {(U • V)Vp' + [(Vp') • V] U} 

+visc, (3) 

where vise gathers all the term containing v or v t . This equation is only an approximation, in so far as 
the assumption of " rapidly varying scales" becomes less and less valid as the resolved scales and subfilter 
velocities are becoming closer and closer in scale space. However, it seems to capture the dominant physics 
of the evolution of the vector 1, as will be later shown. Further, it is tempting to simplify the viscous 
terms of the second equation of J3J to try and get a closed equation for 1. Indeed, these terms involve an 
a priori rather arbitrary turbulent viscosity and one could redefine it so that the viscous terms are simply 
lumped into a term V(^ + v t )V\. Finally, we note that the terms involving the subfilter pressure depend 
on boundary conditions and on subfilter velocities, so that they can be thought to vary over a fast time 
scale, contrarily to the f term which varies over a slow time scale. We therefore choose to collect all the 
term involving the pressure and f into a noise term , £o + £, such that £ is a Gaussian-centered noise and 

£o = -(U-V)f-(f-V)U, 

<£>=0, 

<U{t,x)^{t',x') >=T ij (t-t',x-x'), (4) 

where is the noise correlation function, to be specified. Note that since £ comes from a pressure 
contribution, it only affects the non-solcinodal part of the turbulent force. In the sequel, since we work with 
periodic boundary conditions, we can apply a simple projection procedure to work only with soleinodal 
fields, thereby discarding the £ term. We therefore leave the study of the influence of this term for future 
work, involving more realistic boundary conditions such as flow near walls. 

Collecting all the results, we therefore obtain the following RDT based-model for the turbulent force 
1 as: 

d t l = -(U • V)l' - (1' ■ V)U + VO + i/ t )VI + £o + £ (5) 

It takes the form of a generalized Langevin equation, with friction made of viscosity and rapid distortion 
by resolved scales, and with stochastic forcing, with mean value generated through the energy cascade, 
and whose correlations are physically imposed by boundary conditions via pressure terms. 

Before implementing this Langevin equation into a LES model, we first validate it through dynamical 
a priori tests and seek optimal performances by tuning of v t . Our LES-Langevin model will follow after 
this validation step. 

2.2 Dynamical a priori testing 
2.2.1 The numerical procedure 

The numerical simulations described in the present paper have been performed with a spectral code over a 
cubic domain. This choice is mainly motivated by its performances in terms of accuracy and its versatility 
regarding any variant of our model. The aliasing is removed by keeping only the 2/3 largest modes in 
each direction. The time integration is performed with a second order Adams Bashforth scheme. The 
simulations with both the RDT and the LES-Langevin models are performed with a sharp cut-off filter 
at k = k c . The filter may not be optimum in terms of efficiency, but it has the interesting properties to 
clearly separate scales and to commute with derivative operators. The influence of the filter would need 
to be conducted for a complete validation. The validation tests are performed on decaying and forced 
isotropic homogeneous turbulence. The simulation of decaying turbulence is initialized with a Gaussian 
velocity field with a given energy spectra E(k, t a ) oc k 4 exp(—k 2 /8). For the forced simulation, the forcing 
F is defined by F(k, t)dt — 7(fc,t)u(k,i), with 7 chosen so that the difference of the energy spectra 
before and after the forcing (i.e. the energy injection rate 3) is constant in time (3 = 1.). The forcing is 
concentrated at the lowest wavenumbers k such as < k < 1. The initial time of the forced simulation is 
chosen so that statistical properties are stabilized. The simulations are integrated over approximately 6.3 
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molecular viscosity of v — 0.0014 is used for all simulations. The resulting Taylor Reynolds number is 
approximately constant and equal to R\ = 200 for the forced DNS and varies from R\ — 260 to R\ = 26 
in the decaying case. Both DNS are performed with 341 3 effective wavenumbers after desaliasing (512 3 
grid points). The space resolution allows a computation of scales down to 0.5 and 0.8 Kolmogorov scale for 
the decaying and the forced DNS respectively. The simulation with LES-Langevin model are performed 
with k c — 21 (42 3 effective wavenumbers are used for the resolved scales). 

2.2.2 Validation of the model 

The validation of the model is conducted by dynamical tests performed with the DNS of decaying 3D 
turbulence. It is performed by comparing a full DNS and a simulation at the same resolution, in which 
the turbulent force 1 is replaced by the solution of ©. This way, we may explore the validity of the 
approximations we made in its derivation. We focus here on the energy spectra, so as to capture possible 
deficiency of our model regarding energy transfer between scales. The test is performed for 1 < t < 2 (all 
the simulations are initialized with the same field at t=l). By that time, the kinetic energy is divided by 
a factor of approximately 2. The energy spectra at time t — 2 are shown in fig. ^ 




Figure 1: Energy spectra of the RDT Models © and the equivalent DNS at t = 2 of simulations of decaying 
turbulence (the RDT simulation is initialized with the velocity field of the DNS at t=l) 

The model (RDT) and the DNS agree at the largest scales (k < 3) but they significantly differ for 
smaller scales (the bump at k=21 is due to the coupling between the equation for resolved scales and 
the equation for subfilter scales running in parallel). This can be explained through the analysis of the 
evolution of 1 (fig. |5J). 

One sees that the RDT model leads to a constant increase of the smallest modes of 1 in time. After a 
given period of time, the contribution of these unrealistic resolved scales of 1 influence the model of the 
velocity field resolved scales. 

An explanation of this feature can be found following a recent study of the RDT model (0) |T2*] 
showing that the process of small-scale stretching by random large scales is akin to a dynamo process, 
with exponential increase of the small-scale energy. A way to stabilize the system is to include a friction 
term in the RDT equation, leading to a stationary energy spectrum with index depending on the friction 
time Tf. A Kolmogorov fc -5 / 3 spectrum is obtained for r/ = 27/22 fl, where fl = ((SijSij)^ 1 ^ 2 ) is a 
typical stretching rate based on spatial average of the large-scale velocity stress tensor Sij. Using eq. 

one sees that such a friction term generates an equivalent friction term in the equation for the Lamb 
vector. This remark motivates the introduction of a stabilizing friction term — 1/r/ in the equation for 1 
to try and stabilize the coupled system. Indeed one observes a significant improvement with respect to 
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Figure 2: Comparison at two different times of the spectral density of l 2 modeled by the integration of 1 (eq. 
[SJ and the same quantity directly computed from the resolved and sub-grid scales (1 = (U- V)u' + (u'- V)U). 

the original RDT model. We therefore adopt this procedure as our starting RDT model, from which we 
now build our Langevin-LES model. 

2.3 Derivation of the model 

The derivation of the LES-Langevin model of turbulence proceeds in two steps. In the first one, we 
replace the term (u' • V)u' into a turbulent viscous term /j, t AU, acting only at large scales, in the spirit of 
standard deterministic " eddy- viscosity" models. In a second step, we derive a suitable Langevin equation 
of the turbulent force 1 through a decimation of the number of degrees of freedom corresponding to scales 
beyond the aliasing limit. For this, we introduce a strong hyperviscosity to damp all components beyond 
a given cut-off wavenumber k m . There is a priori complete freedom for the choice of k rn . Here, we note 
that the cascade-driven forcing f has components only up to k = 2k c . Therefore, any component of 
1 beyond 3fc c will only be generated through secondary processes (stretching) rather than through the 
forcing. We may then hope that any k m between 2k c and 3fc c provides the dominant contribution to the 
stochastic term 1. We found that k m = 2k c is in fact sufficient to capture this dominant contribution. 
Our Langevin-LES (LRDT) model is therefore finally given by: 



d t XJ + (U ■ V)U + 1 = -VP + V(z/ + /i t )VU; [0 < k < k c ] 
d t l = -(l/T f )-(U-V)l'-(l'-V)U 



+V{u + u t )Vl + Co + f ; [0 < k < 2k c ) (6) 
£o = -(U-V)f-(f -V)U, 

< £ >= 0, < Ci(t,x)Cj(t',x') >= F i: j(t-t',x- x'), 



where r/ = 27 /22((SijSij) 1 ^ 2 ) , v t and [it will be specified later and fi = dj (UiUj — UJJj). Looking at eq. 
10, one recognizes a LES model where the backscatter coming from resolved scales-subfilter interaction 
is parameterized through a noise. The latter obeys a generalized Langevin equation, with friction made 
of viscosity and rapid distortion by resolved scales, and with stochastic forcing £o + generated through 
the energy cascade and pressure processes. 



5 



The RDT-based model we described is fully dynamics, and a priori adapted to any type of geometry, 
with possible anisotropies. In certain very simple cases, however, this model may support additional 
approximations, enabling a reduction in the computation cost. We present two examples below. 

2.4.1 Quasi-linear approximation (LQL) 

In the spirit of quasi-linear approximation, we may lump the transport and stretching by resolved scales 
with the viscous and friction term and simply replace them by a total friction term — 1/r, where r is a 
typical time scale to be chosen later. 

Note that this simplification somehow relies on the hypothesis that the advection and stretching of 1 
by the resolved scales proceeds in an isotropic manner. We therefore do not expect this simple procedure 
to be fully valid in cases with e.g. stratification, or rotation, or walls. In the present case, it appears to 
give very good results, provided the time scale is chosen appropriately. This model therefore looks like: 



dtU + (U • V)U + 1 = -VP + V(f + Mt)VU; [0 < k < k c ] 

8 t l = -(l/r) + Co + £ [0 < k < 2k c ] (7) 

2.4.2 The over-damped approximation (LQLOD) 

In situation where the time scale r is very small, the damping in the Langevin equation is very large, and 
there is almost instantaneous adjustment of 1 to the forcing. This situation typically arises in locations 
where the resolved gradients are very large, i.e;. at the border of coherent resolved eddies. It is therefore 
interesting to run the model in this over-damped limit, to see whether these kinds of event dominate the 
overall dynamics. The resulting model (Langevin Quasi-Linear Over-Damped) looks like: 

d t U + (U • V)U + 1 = -VP + V{v + Mt)VU; [0 < k < k c ] 

1 = r (Co + ; [0 < k < 2k c ] (8) 

Through this approximation, our model, based on the estimation of 1, becomes equivalent to the ADM 
velocity-estimation model of Domaradzki and collaborators IE]- In a recent improvement of this 
model, the estimated small scales are used to integrate the Navier-Stokes equation on the fine grid over 
a given period of time short enough to prevent the pile up of the energy at the finest scales. The ADM 
model can therefore be viewed as a special case of our Langevin model. In the following, we show that 
this approximation perform very well and is indistinguishable from the general model in the isotropic 
homogeneous case. However, it may not be the case in more general situations, such as near walls. We 
leave this for future work. 

2.5 Parameters 

Our model includes three parameters T^, r and fit, that we now discuss. As previously discussed, the 
stochastic correlation comes from pressure contribution, and we found its influence be irrelevant in 
the case with periodic boundary condition. In the present study, we therefore adopt Ty = 0, leaving its 

study for cases involving walls. For the time scale, we set r = t$ + a (SijSij) , so as to allow spatial 
variation of the friction time scale in agreement with local stretching. Alternatively, one may also use a 
global friction time by considering only spatial average of the local stretching rate, so that r = r m = /3r/. 
This last prescription is actually better in term of computational cost. Empirically, we found that both 
prescriptions lead to a good model. In our numerical tests, we chose a = 1 or j3 = 1/2 but we checked 
that the results are not very sensitive to these parameters 

For the turbulent viscosity, we have an a priori rather vast choice, ranging from spectral to Smagorinsky 
models. In the present case, where our numerical scheme uses fast Fourier transform, the optimal choice 
(in term of computational cost) is the spectral model. The simplest choice is to link the turbulent viscosity 
to the kinetic energy at the cut-off wavenumber k c (/i t oc y/ E(k c )/k c ) . The proportionality factor is tuned 
in our model and our numerical tests. Since part of the dissipation is produced by 1, the proportionality 
constant turns out to be reduced from 0.267 for the original spectral model (see j2J) to 0.08 in our case. 
Summarizing, we get the following parameters for our models: 
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Figure 3: Energy spectra of the LRDT model (dash lines) and the equivalent DNS (continuous lines) of a 
decaying homogeneous isotropic turbulence. The spectra are compared for 4 different times: t = (1,3,6,9). 



Tij = (9) 

r = (27/22)((S^-r 1/2 ) + (S«S«)~ 1/2 (10) 
r m = 0.5<(,%,%r 1/2 ) (11) 
IH = 0.08(E(k c )/k c ) 1/2 (12) 

Models using the second prescription for the time scale will be called Langevin Quasi Linear Averaged 
(LQLA). They are optimal in term of computational cost. 



3 Model performances 

3.1 Comparison of the Langevin models 

As a preliminary assessment of the interest of our model, we choose to test it in a simple configuration 
where we have both the code and the competence to implement it over a short time scale. This case corre- 
sponds to homogeneous, isotropic turbulence with periodic boundary conditions. It is clearly a restricted 
problem, where backscatter has proved to be secondary and where eddy-viscosity based approaches per- 
form very satisfactorily. In some sense, this ideal case is therefore a challenge to our model, for it is not 
clear whether there is room for further improvement with respect to eddy-viscosity description. In the 
following, we demonstrate a slight but discernible positive effect of our strategy. 

The test is done both on decaying turbulence and forced turbulence in order to check the model abil- 
ity to accurately reproduce the energy cascade from resolved to unresolved scales. The first comparison 
between the LRDT model (eq. |5J) and the DNS is conducted over the energy spectra in fig. [5] and the 
decay of the total kinetic energy in fig. 0] In this comparison, the resolved scale turbulent viscosity /x t 
is tuned to give the best result in terms of energy spectra and decay of kinetic energy. The decimation 
of the small-scale of 1 (so as to keep only 2k c modes in each direction) is done through an hyper- viscosity 
(y t = 0.002 k 2 y / E(k c )/k c ) super-seeding the turbulent viscosity of the original RDT model. The energy 
spectra are in good agreement with the DNS for t < 1 and slightly differ for the smallest resolved scales 
at later time. We are confident that a more elaborate tuning of the viscosity and friction term could have 
improved the model over later time, but such accuracy was not needed in the sequel. 
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Figure 4: Comparison of the decay of kinetic energy between the DNS of decaying turbulence, the LRDT 
model (with k c = 21, k m = 2k c ) and the Smagorinsky model with k c = 42. The time t = 10 correspond to a 
non-dimensional time of 6.3 

The same numerical tests have been performed for the approximate models, LQL and LQLOD, (figs. 
CDandinj), with the same definition of the time scale r (eq. JUJ , dJ ) ■ The results are slightly less accurate 
than the original LRDT model, but the two approximated models represent a good compromise between 
accuracy and CPU cost. As discussed in the previous section, an alternative model to LQL can also 
be formulate using an averaged value r m of the time scale t. A comparison of averaged energy spectra 
between LQL and LQLA models and DNS is done fig. 0for a simulation of forced homogeneous isotropic 
turbulence. In the simulation with LQLA model, the averaged value of the time scale r m is of the order 
of 0.1 turnover time. The two LES simulations are very similar and are both in good agreement with the 
DNS. 

3.2 Energy transfer 

One of the ability of the model compared to eddy viscosity model is to allow energy transfer from subfilter 
scales to resolved scales through the model of the cross stress tensor by the Langevin equation of 1. The 
amount of backscatter can be quantify by: 



' res, sub 



00=/ |u(p)U(q) ( 5(k-p-q)d 3 pd 3 q (13) 

The integration of T reSjSU b(\i} over a shell leads to the energy transfer as a function of the k. This quantity 
was investigated for the LQLA model in the case of forced isotropic turbulence (see fig. |SJ. The intensity 
of the energy transfer through the Lamb vector T^®^ b (k) is compared with the energy transfer from 
resolved scales T^Q^ea s ) anc ^ the energy dissipation introduced by the model of the Reynolds stress tensor 
with a turbulent viscosity T^~^ b (k) . Each component of the energy transfer is also compared with the 
equivalent component computed from the DNS. The statistics have been conducted over 4 turnover times 
after stabilization of the statistics. The corresponding average energy spectra of the two simulations are 
shown in Fig. El Since the two solutions used for the statistics differ in time, it is difficult to compare 
exactly the transfer between the model and the DNS. However one can compare the sign and the order 
of magnitude of each component. The transfer coming from the Lamb vector exhibits a small amount 
of backscatter, which is not present in the DNS. However the transfer is very comparable at higher 
wavenumbers. The turbulent viscosity /i t used to model the Reynolds stress tensor leads to a dissipation 
of energy higher than the equivalent energy transfer from the DNS. The model of this term could probably 
be improved. 
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Figure 5: Energy spectra of the LQL model (dash lines) and the equivalent DNS (continuous lines) of a 
decaying homogeneous isotropic turbulence. The spectra are compared for 4 different times: t = (1,3,6,9). 




Figure 6: Energy spectra of the LQLOD model (dash lines) and the equivalent DNS (continuous lines) of a 
decaying homogeneous isotropic turbulence. The spectra are compared for 4 different times: t = (1,3,6,9). 
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Figure 7: Averaged energy spectra for the simulation of forced isotropic homogeneous turbulence. The 
spectra are averaged over 4 turnover times for the DNS and 40 turnover times for the three LES. The DNS 
was performed with 384 3 Fourier modes. 




Figure 8: Comparison of the energy tranfer between the DNS and the LQLA model for forced isotropic 
turbulence. The energy tranfer is split into the 3 components involving resolved and subgrid scales. T^^ b 

is the energy transfer from the Lamb vector (see eq. IT3]) and T^2^ ub is the disipation introduced by the 
turbulent viscosity fit- 
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Figure 9: Comparison of the energy spectra of our LQLA model with the DNS, the Smagorinsky model 
and the spectral model (averaged over 4 turnover times for DNS and 40 turnover times for the 3 models) 



3.3 Comparison with other LES models 

Within the present framework (isotropic, homogeneous turbulence) it therefore seems that the LQLA 
model represents the best compromise in terms of accuracy and CPU cost. We now explore its perfor- 
mances with respect to other classical LES models in this field. 

The simpler LES models for homogeneous, isotropic turbulence are based on eddy-viscosity concept, 
namely the Smagorinsky model [3^ and the spectral model |23. The later has only been tested for forced 
isotropic turbulence. The former corresponds to a Smagorinsky model (without dynamical procedure) in 
which the constant was tuned to get the best results in terms of energy decay. It was run on the finest 
grid of our model (the grid used for the integration of 1) . The comparison is done using several statistical 
diagnostic, tracing different properties of the model. 

3.3.1 Energy spectra 

A comparison of the energy spectra for these two models and our LQLA model is shown in fig. 1101 for 
the decaying case and in fig. [5] for the forced case. In the decaying case, our Langevin model performed 
at least as well as the Smagorinsky model in terms of energy decay. In the forced case, the Smagorinsky 
model clearly overestimate the energy dissipation at the smallest resolved scales, resulting in an energy 
spectrum much steeper than the theoretical —5/3 slopes. Similar result is found for the spectral model, 
while the LQLA model compares much better with the DNS at all scales. Indeed, our model seems to 
model properly the energy transfer at the largest scales and the energy dissipation near the cut-off. 

3.3.2 Skewness and flatness 

We also compared the skewness and the flatness of the LES models to the same statistics for the corre- 
sponding scales of the DNS in the case of forced turbulence. The results are displayed in tabled The 
LQLA model gives a correct sign of the skewness, with a value three times too large, and a correct esti- 
mate of the flatness. By comparison, the Smagorinsky model predicts a wrong sign for the skewness and 
underestimates slightly the flatness. The spectral model gives the best estimate of the skewness (with 
correct sign) and underestimates the flatness of the distribution. 
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Figure 10: Energy spectra of the LQLA model (dash lines) and the equivalent DNS (continuous lines) of a 
decaying homogeneous isotropic turbulence. The spectra are compared for 4 different times: t = (1,3,6,9). 

Table 1: Skewness and Flatness of the resolved scales of the LQLA model, the spectral model, the Smagorin- 
sky model and the corresponding DNS in the case of forced turbulence. 



model 


skewness 


flatness 


DNS 


-0.003314 


0.3141 


LQLA 


-0.011198 


0.3107 


Spectral 


-0.004328 


0.3094 


Smagorinsky 


0.013744 


0.3078 



3.3.3 Probability distribution function 

A more refined statistical comparison can be performed using the Probability Density Function of the 
velocity increment 5v r = (v(x + r) +v(x)) (figs. IHIand U^! . Because of the coarse resolution of the LES 
simulations and the small number of statistics for the equivalent scales of DNS, the differences of PDFs 
between the models and the DNS are difficult to observe. However, the statistics with the LQLA model 
seems to be in a closer agreement with the DNS for both longitudinal and transverse velocity increment 
than the two other dissipativc models. 

3.3.4 Statistics of invariants 

A different statistical test can also be performed on the tensor of velocity gradient tensor Aij = duljdxj. 
We choose to focus on the two following invariants: 

Q ~^Ai n iA rn i (14) 

R — ~-^Ai rn A rn ]~Aki 1 (15) 
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Figure 11: Comparison of the normalized PDF of longitudinal velocity increments (6v r = v(x + r) — v(x) 
with |r| = ^ and r x v = 0) for our LQLA model, the Smagorinsky model, the spectral model and the 
DNS. The PDF are averaged over 4 turnover times for DNS and 40 turnover times for the LES models. 




-5 
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Figure 12: Comparison of the normalized PDF of transverse velocity increments (6v r = v(x + r) — v(x) 
with |r| = |j| and r • v = 0) for our LQLA model, the Smagorinsky model, the spectral model and the DNS. 
The PDF are averaged over 4 turnover times for DNS and 40 turnover times for the LES models. 
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Figure 13: Comparison of the joint PDF of R* and Q* (eqs. ITT)) for the resolved scales of the LQLA model 
(upper right), the Smagorinsky model (lower right), the spectral model (lower left) and the corresponding 
DNS (upper left) in the case of forced turbulence. 



which have been used to discriminate between different turbulent models . For this, we normalize the 
two invariants with a time scale based on the strain-rate tensor. The two invariants become: 



(SijSij) 



where s.y = Sij — (Sij) and (•) is a sample averaging. The figure 1X31 shows the joint PDF of Q* and R* 
for the resolved scale of the LQLA model, the Smagorinsky model, the spectral model, and the DNS. One 
sees that the LQA model and the spectral model reproduce better the joint PDF than the Smagorinsky 
model, which tends to squeeze the level lines towards R — 0, a feature already noted in |38|. 



(16) 
(17) 
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Figure 14: Comparison of the PDF of s* (eq. I18J) for the resolved scales of the LQLA model, the Smagorinsky 
model and the spectral model and the corresponding DNS in the case of forced turbulence. 



3.3.5 Strain statistics 

Another important characteristic of the strain rate tensor S is the statistics of its eigenvalues. To quantify 
it, Lund and Rogers [U] proposed the following parameter: 

»• = -^j^m < 18 > 

where a, j3 and 7 are the eigenvalues of S such as a < (3 < 7. The term s* measure the shape of the 
deformation caused by the strain rate tensor, with s* = —1,0,1 corresponding respectively to worms, 
shear and pancakes. The probability density function of s* for the resolved scales of the LQLA model as 
well as for the Smagorinsky and the spectral models arc computed in the case of forced turbulence. The 
comparison is shown on figure El The comparison shows very good agreement for both the LQLA model 
and the spectral model with respect to the DNS, while the Smagorinsky model shows poorer agreement. 



3.3.6 Computational cost 

The performance of each model can also be analyzed with respect to its computational cost. However, 
we are not able to give precise estimations of the performance because the implementation of each model 
was not optimized in terms of CPU cost. Moreover, the cost is highly related to the type of discretization 
and the numerical scheme. But, looking at the equations, we can expect an optimal computational cost 
for our model to be of the order of a DNS on a grid corresponding to the grid used for the integration 
of 1. It as been shown that an integration with wavenumbers up to 2k c is enough for the equation of 1. 
With the over-damped version of the model, there is no integration for 1 and we could expect a gain of a 
factor 2 relative to the time step. 



4 Discussion 

In this work, we present a new strategy to model subfilter scale in a Large-Eddy-Simulation. Our phi- 
losophy derives from the Brownian motion, where the extra degrees of freedom are modeled through the 
combination of a deterministic friction and stochastic forcing. Our model therefore includes both a de- 
terministic eddy-viscosity and a stochastic turbulent force, solution of a generalized Langevin equation. 
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motions, and stochastic energy backscatter from subfilter to resolved scale. Several model have been 
proposed in the past, to reproduce these two effects. At this point, it is interesting to review briefly the 
most popular ones, so as to understand their limitations and differences with respect to our model. A first 
example is the similarity model, which does not assume the alignment of the turbulent stress tensor with 
the strain rate tensor. In the original model of Bardina £Q, the model of the Reynolds stress is achieved 
through approximation of the unfiltered velocity by the filtered velocity. This model can be seen as the 
zero-level of more general deconvolution models that rebuilds the unfiltered velocity (see for a dis- 
cussion on de-filtering in LES). This procedure however is also subject to serious limitations, since a LES 
with an exact deconvolution of the velocity field would lead to a simulation of the Navier-Stokes equation 
on a coarse grid without any model but with a sharp filter. Practically, the resulting lack of dissipation is 
provided through the physical or the numerical approximations of the deconvolution . Other concepts 
arose from a detailed analysis of the energy transfer between resolved and subfilter scales. This leads to 
the conclusion that the scales smaller than half the smallest resolved scales have only few impact on the 
direct energy transfer to the largest scales 6, 40, 4 . Following this observation, a new method to capture 
the backscatter has been proposed through direct estimation or a more accurate integration of the largest 
subfilter-scales velocities. The simulations with these models can be accurate but expensive LES or cheap 
approximated DNS. The dynamics multilevel model |l()l is an example of cheap DNS as it requires 

the same resolution as the corresponding DNS. Our 2D version of the RDT model using a Lagrangian 
evolution of small-scale wave-packets is more flexible since the accuracy of the model is function of the 
number of modes used for the subfilter scales. Better performances in terms of computational cost can be 
obtained in models where only the largest unresolved scales are estimated. An example is provided by the 
ADM model of Domaradzki 00 IE], where subfilter scales are estimated on a finer grid. The estimation 
of subfilter scales is performed after a first step of deconvolution. In a recent improvement of this model 
the estimated small scales are used to integrate the Navier Stokes equation on the fine grid over a given 
period of time short enough to prevent the pile up of the energy at the finest scales. Finally, one may note 
that attempt of direct inclusion of the backscatter have been done in the past through empirical addition 
of random numbers, with well chosen spectrum |2()l I25| . These models led to noticeable improvement 
with respect to eddy-viscosity type of approach for boundary layer or plane shear mixing layer. 

Most of these models have been validated and compared in several flow configurations. The most 
intensive comparison have been performed for isotropic turbulence where the DNS can reach higher 
Reynolds number. For instance, Fureby et al ^Sj made a comparative study of subgrid scale model of 
eddy-viscosity type, models of similarity type and Miles model. For moderate Reynolds number, they 
concluded that the difference between LES with different models are small but not insignificant. The 
characteristics of more specific type of models like eddy-viscosity models or similarity models have been 
studied intensively by many authors f |2"51 12"?1 12~2"] 

From this review, it is clear that our originality is to seek for a systematic derivation of the stochastic 
effect, rooted in the Navier-Stokes dynamics. By this mean, we hope both to avoid uncontrolled, empirical 
modeling, and respect of all the symmetries of the original Navier-Stokes equation (a constraint not always 
easy to respect, see |29p. Our model may also easily be generalized to other more complicated systems, 
like e.g. rotating, stratified or inhomogeneous turbulence. Specifically, we separate the subfilter scale 
contribution into a term that correlates well to resolved strain tensors, and a term susceptible to the 
backscatter. The first term is modeled through a traditional eddy-viscosity, while the second is modeled 
through a noise, obeying a generalized Langevin equation. This equation is not postulated, but derived 
from the Navier-Stokes equation, through hypothesis akin to rapid distortion theory. The stochastic retro- 
action couples the resolved scales and the noise through the energy cascade, ensuring strong non-linearity 
of our model and non-trivial equilibrium solutions. A corner stone of our model is the estimation of 1, 
rather than the velocity, like in ADM model. Apart from this difference, there is some similarity between 
our model and the ADM model, in the simplest approximations we considered (quasi-linear and over- 
damped approximation). However, our model derives from the original Navier-Stokes equation through a 
well-defined and systematic procedure. 

We have proposed a preliminary test of our model in isotropic, homogeneous situation, with periodic 
boundary conditions. In this case where the backscatter is secondary, we have shown that our model 
still perform slightly better than traditional LES approach used in this situation, based on eddy-viscosity 
concept. It would now be interesting to test our model in more realistic situations, with boundaries, 
where the stochastic modeling will become an essential ingredient. Given that our model reduces to the 
ADM model in a well-defined limit, we are confident that it can perform as well as this last model. The 
interesting question is whether our dynamical procedure increases its performance, and overcome certain 
limitations of the ADM model. This will be the subject of a future work. 
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